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The 3-level leapfrog time integration algorithm is an at- 
tractive choice for numerical relativity simulations since it is 
time-symmetric and avoids non-physical damping. In Newto- 
nian problems without velocity dependent forces, this method 
enjoys the advantage of long term stability. However, for more 
general differential equations, whether ordinary or partial, de- 
layed onset numerical instabilities can arise and destroy the 
solution. A known cure for such instabilities appears to have 
been overlooked in many application areas. We give an im- 
proved cure ("deloused leapfrog") that both reduces memory 
demands (important for 3 + 1 dimensional wave equations) 
and allows for the use of adaptive timesteps without a loss 
in accuracy. We show both that the instability arises and 
that the cure we propose works in highly relativistic prob- 
lems such as tightly bound geodesies, spatially homogeneous 
spacetimes, and strong gravitational waves. In the gravita- 
tional wave test case (polarized waves in a Gowdy spacetime) 
the deloused leapfrog method was five to eight times less CPU 
costly at various accuracies than the implicit Crank-Nicholson 
method, which is not subject to this instability. 
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I. INTRODUCTION 

Numerical relativity comprises the dynamical solution 
of the Einstein equations on a computer, allowing the 
construction of spacetimes that cannot be studied by 
purely analytic methods. A major application of numer- 
ical relativity is the modeling of astrophysical sources of 
gravitational radiation such as binary black hole § or 
neutron star inspiral Q , and nonspherical stellar collapse 
H . The continued development of gravitational wave de- 
tectors, with the expectation that ground-based interfer- 
ometers such as LIGO (§, VIRGO § and GEO600 g 
will begin taking data in a few years, gives these studies 
a high priority. Numerical relativity is also important for 
studying the dynamics of pure gravitational waves , in- 
homogeneous cosm olog ies ||, the behavior of cosmolog- 
ical singularities p|Jl0||, and critical behavior in general 



relativity Ipif. 

All of these endeavors require accurate numerical al- 
gorithms to correctly model the physics of curved space- 
time. Simulations in three spatial dimensions plus time 
are expensive in terms of both CPU usage and memory 
requirements, and thus demand numerical methods that 
are efficient in both these regards. Memory limits, how- 
ever, are less elastic in the short term than CPU time 
constraints. Thus a three-level second order algorithm 
may be more appropriate than a faster, high order algo- 
rithm which can only be implemented on smaller prob- 
lems. Also, modeling the inspiral of binary black holes 
or neutron stars requires evolving the system for many 
orbital periods, so that numerical algorithms with long 
term stability and freedom from unphysical damping are 
essential. 

Leapfrog methods are often used for the time inte- 
gration of equations in numerical relativity and other 
branches of computational physics. The 3-level leapfrog 
method has the important property of being symplectic. 
In the context of a Hamiltonian system for which the 
differential equation has a symplectic structure (conju- 
gate pairing of coordinates and momenta), this means 
that the difference equations also have such a structure 
and the integration step in the difference equations is a 
canonical transformation. With a symplectic integrator, 
all the Lagrangian integral invariants, including phase 
space volume, are exactly conserved by the integration 
scheme. Since the leapfrog method is time symmetric 
and maintains good conservation of physically conserved 
quantities |U^-|l4|], it has a well-deserved reputation in 
the context of Newtonian mechanics. Unfortunately this 
reputation is generally not merited when velocity depen- 
dent forces are met. In the integration of systems with 
such forces, this scheme is well-known to be suscepti- 
ble to numerical instability (e.g., fil^-|l9[|, and references 
therein), even under conditions where local linearization 
analysis anticipates stability. This instability occurs in 
the integration of both ordinary and partial differential 
equations and, in the case of partial differential equa- 
tions, is independent of the mesh size used for the spatial 
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discretization |17| ]. 

An understanding of the origin of this instability 
was given by Sanz-Serna [ p"7| who pointed out that 
the leapfrog scheme approximates not merely the in- 
tended differential equation system but a larger "aug- 
mented" system containing additional, nonphysical, par- 
asitic modes. Since the leapfrog method is symplectic 
as applied to the augmented system |^0| , the advantages 
of symplectic methods (see |Q) would be attained to 
the extent that the parasitic modes remain zero numer- 
ically, as they do in an exact solution of the augmented 
system pl| . Aoyagi and Abe p8p^ ] identified the diag- 
nostic symptom of this instability as a sawtooth oscilla- 
tion or alternation of values between odd and even steps 
of the integration, and supplied a cure — a Runge-Kutta 
smoother to suppress this oscillation. Subsequent work 
||l9f shows these phenomena in ordinary differential equa- 
tions, where the delayed onset of this instability is clearly 
apparent. 

We have studied the use of the 3-level leapfrog method 
in numerical relativity. In this work (see also we 
extend the ideas of Aoyagi and collaborators for remov- 
ing the unstable parasitic modes, yielding an algorithm 
that reduces the number of time levels of data that must 
be stored by the code, allows the timestep to be changed, 
and thus is better suited to long-term integration of large 
scale numerical systems. Although we concentrate on 
the ADM 3+1 formalism for numerical relativity, our 
methods are quite general and thereby applicable to a 
wide range of problems in computational physics. Sec- 
tion [n| describes the source of the 3-level leapfrog instabil- 
ity. Section III presents our algorithm, dubbed "deloused 
leapfrog," for removing instabiliti es th at may arise in 3- 
level leapfrog integrations. In Sec. IV, we summarize the 
3 + 1 formalism of numerical relativity and in Sec. |v| 
introduce three simple classes of Gowdy T 3 spacetimes. 
In Sec. fv|, we use these models as relativistic testbeds 
for the deloused leapfrog technique and evaluate the nu- 
merical efficiency of deloused leapfrog by comparing its 
cost-effectiveness with those of the staggered leapfrog and 
Crank-Nicholson techniques. A summary and discussion 
of the results of the stability and efficienc y tes ts of the 



deloused leapfrog technique is given in Sec. VII 



II. THE 3-LEVEL LEAPFROG INSTABILITY 

We begin by considering the system of differential 
equations 



dz 



F(z,i). 



(1) 



Although we use a system of ordinary differential equa- 
tions [Eqs. (|l|)] to describe the leapfrog instability and 
its cure in this section and the next, the techniques we 
outline apply equally well to systems of partial differen- 
tial equations in which the time integration is carried out 



using the leapfrog method; see Sec. VI below. The 3-level 
leapfrog discretization of this system is 



2F(i n+1 ,t n+1 )At, 



(2) 



where we assume a constant timestep At between time 
levels n and n+1. The distinction between z in the differ- 
ential equation and z in the difference equation is a warn- 
ing that the relationship is not as straighforward as first 
appears. The discretization given in Eqs. (|) has 0(At 2 ) 
accuracy and is a 3-level method, in that knowledge of 
data on time levels n and n + 1 is needed to compute the 
result on time level n + 2. The 3-level leapfrog algorithm 
[Eqs. (||)] is symplectic and time reversible, which means 
that it provides a Hamiltonian (damping free) model of 
an underlying Hamiltonian differential equation system 

us 

Equations (pi) arise in many physical applications in 
which z comprises both position and velocity data z = 
(r , v) . For example, the motion of a particle of mass /i 
under the action of a force / = fid is given by the set of 
equations 



dr 
~dl 
dv 
~dt 



(3) 



where r is the position vector of the particle. The 3-level 
leapfrog discretization of Eqs. ([|) gives 



r n+1 = r + 2v n At 
v n+1 =v n ~ 1 + 2a n At 



(4) 



For a particle moving in a Newtonian gravitational field 
a n = a n (r n ) and the system given by Eqs. (||) is stable 
and thus suitable for long time integrations. However, if 
there are so-called "velocity dependent forces" in which 
a" = a n (f v n ), such as arise for a particle moving un- 
der the influence of a magnetic field or a general relativis- 
tic gravitational field, the 3-level leapfrog scheme [Eqs. 
(^)] can become unstable. As we shall explain below, 
nonphysical parasitic modes can arise during the time in- 
tegration and eventually destroy the numerical solution. 

Notice that the leapfrog algorithm [Eqs. (||)] gives the 
value of z™ +2 at, say, the even time level n + 2 in terms 
of the value z™ at the even level n and the source term 
F(z n+1 , t n+1 ) at the odd level n + 1. It also requires the 
specification of initial data at two time levels, z° and z 1 , 
which is twice as much as the original first-order differ- 
ential equation system requires. This doubling of initial 
condition specifications is the clue to the fact that this 
numerical algorithm using z has twice as many degrees 
of freedom as does the physical system where states are 
specified by z. The algorithm can be expressed, with a 
change of notation, by writing the solution at the even 
timesteps as z 2n = x 2 ™ and at the odd timesteps as 

~2n+l = y 2n+l W j th ^ Eqg (j| become 
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,2n+2 



Jin 



x— =x-" + 2F(y^ n+1 ,^" +1 )Af 



,2n+3 



,2n+l 



2F(x 



2n+2 ^.2n+2 



(5) 



Sanz-Serna [|17| (concisely summarized in |2j|) notes that 
Eqs. dl) can be considered to be a consistent single-step 
discretization of a larger system of equations for the even 
solutions x and the odd solutions y, 



(G) 



with timestep 2Ai. Equations (g) are known as the aug- 
mented system ]l7| , p3| ]. 

If z is a solution of Eqs. (|l|), it gives a solution of 
Eqs. (^), the augmented system, as x = z and y = z. 
In general, however, other (unphysical) solutions will be 
possible. Since these unphysical, parasitic solutions can 
arise as valid solutions to the augmented difference sys- 
tem of Eqs. (|5|), the 3-level leapfrog integrator alone can- 
not distinguish between them and the physical solutions. 

Since the physical solutions of Eqs. (g) have x = y it 
is natural to define 

x = z + w , y = z w , (7) 

so that the w measure the parasitic deviations from the 
desired physical solutions z. One can rewrite Eqs. (Q) 
in finite difference form, using the definitions preceding 
Eqs. (El), as 



-l)"w r ' 



(8) 



The solution z™ to the leapfrog difference system [Eqs. 
(Q)] thus contains both physical z" and parasitic modes 
w™. With this notation, the augmented differential equa- 
tions can be written as 1119(1 



dz 

2— = F(z - w, t) + F(z + w, t) 
dt 

dw 

2 = F(z - w, t) - F(z + w, *) 

dt 



(9) 



We note that when these are expanded in powers of w 
one has 



^=F(z,i) + 0(w 2 ) , 
^ = -DF(z,t) • w + C(w 3 ) 



(10) 



where DF is the matrix of partial derivatives dFi/dzj. 

The second of Eqs. ([n]), ignoring the cubic and higher 
terms, is a linear equation (in w) which from the time 
dependence in DF(z(t), t) easily gives rise to paramet- 
ric amplification p4[ leading to growth of the parasitic 
modes w. Such parametric amplification was clearly di- 
agnosed in a useful example |18[ ). A still stronger ar- 
gument for the growth of parasitic modes (without the 



above linear perturbation assumptions) had been given 
earlier by Sanz-Serna 0j2^] based on his proof that the 
leapfrog scheme preserves volume in the augmented state 
space even when the original system is not Hamiltonian. 
Sanz-Serna has an elementary example, dz/dt = z 2 , 
z(0) = —1, with a two dimensional zw augmented phase 
space (Fig. I in |rij]); this shows all the qualitative fea- 
tures of the delayed onset instability typical of numerical 
catastrophies that generally result from using leapfrog 
in nonlinear systems with velocity dependent forces. In 
this example, every solution diverges [z — > oo) if at any 
point w ^ 0, although z —> as t — * oo for the physi- 
cal w — solution. Numerical experiments with leapfrog 
follow the divergent solutions of this analytical example 
since w ^ at some point arises either from imperfect 
initial conditions, roundoff error, or discretization error. 
Hence, nonlinear interaction or parametric amplification 
of linear interaction between the physical and parasitic 
modes of the solution to the leapfrog difference system 
[Eqs. (||)] can cause the parasitic mode to grow to the 
point where it destroys the numerical solution. 

However, if the augmented difference system [Eqs. (Q) 
or (||)] decouples, no such instability occurs. To illustrate 
what it means for these equations to decouple, consider 
the augmented difference system for Eqs. (|), in the ab- 
sence of velocity dependent forces. In that case one can 
let x = (r, u) and y = {q, v) to find that the (even r, 
odd v) system does not couple to the (odd q, even u) 
system. Thus the augmented difference system consists 
of two interlaced, noninteracting copies (even r, odd v) 
and (odd q, even u) of the physical system. Each of these 
two systems is of the Newtonian form where the leapfrog 
scheme has shown itself remarkably stable. 

Thus, as mentioned above, in a 3-level leapfrog inte- 
gration of Eqs. (|3|) stablity can generally only be antici- 
pated in the absence of velocity dependent forces. In this 
case, the parasitic mode is generally still present as the 
difference between two interlaced numerical solutions of 
Eqs. (S), but will remain small unless the physical system 
is highly sensitive to small differences in initial conditions 
(i.e., chaotic). Note that it is customary, in the absence 
of velocity dependent forces, to omit the q and u vari- 
ables from the integration scheme, yielding the staggered 
leapfrog scheme (see below). Alternatively, in a code 
based on Eqs. (||), one could in this case solve for two 
independent solutions approximating Eqs. ([!]) based on 
distinct initial conditions for (f,v) and for (q,u). 

A stable 3-level leapfrog integration of a system of 
equations containing "velocity dependent forces" can be 
maintained if the growth of the parasitic mode can be 
controlled. When a constraint w = is adjoined to 
Eqs. (9|) they reduce to the original physical system of 
Eq. (p). This constraint, when imposed initially, is pre- 
served by the differential equation system @ since w = 
gives dvt/dt = 0. But in numerical implementations 
errors will inevitably introduce nonzero w. The cure 
proposed by Aoyagi and Abe [jl^,^2| is to reimpose the 
constraint w = as necessary to suppress the parasitic 
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mode. Their method for doing so is based on the identifi- 
cation of w via its signature even-odd timepstep alterna- 
tion in sign, which is evident in Eq. (||) . As noted in , 
this method is very efficient since, from the power se- 
ries expansions in Eqs. (fL0|), the parasitic modes w need 
merely be suppressed to single precision accuracy to as- 
sure that they do not contaminate the physical solution z 
in double precision. We will return to the cure in Sec. Ill 
below. 

Note that the instability under discussion here is not 
related to the "mesh-drifting" instability inherent to non- 
dissipative leapfrog integrations of partial differential 
equations |2^] . As has been emphasized above, the 3-level 
leapfrog instability can arise in integrations of both ordi- 
nary and partial differential equations and results from 
the temporal, not spatial, discretization of the method. 

To demonstrate this instability, we have used the 
3-level leapfrog method to numerically integrate the 
geodesic equations for a particle moving on a circular 
orbit of radius rg in the Schwarzschild spacetimc. In the 
limit ro 3> M, we recover the usual Newtonian equations 
of motion [Eqs. (|j|)]. In this case, the leapfrog method 
produced a circular orbit that was stable for ten thousand 
orbital periods (before we terminated the run). However, 
as r is decreased, general relativistic effects give rise to 
terms that behave like velocity-dependent forces and the 
integrator fails to maintain a stable evolution. Fig. [j] 
shows the results of the geodesic integration for the case 
ro = 1071/. The particle orbit shown in Fig. |l|a initially 
appears to be stable; eventually the instability manifests 
as the solutions on the even and odd timesteps diverge. 
Fig. [[J) shows the magnitude of the particle's position 
vector as a function of time; again the even and odd so- 
lutions clearly diverge as the parasitic mode grows to de- 
structive levels. Although the parasitic mode is present 
from the beginning of the calculation, it takes a number 
of orbits before it grows to noticeable levels; the instabil- 
ity then grows catastrophically, causing the integration 
to crash after about six orbits. 

Two other integration methods (whose efficiency we 
will compare with that of our deloused leapfrog algorithm 
in Sec. VI), staggered leapfrog and Crank-Nicholson, 
do not suffer from the instability present in the 3-level 
leapfrog algorithm. In both cases this happens because 
neither of these algorithms augments the phase space (or 
state space) of the problem by adding new degrees of 
freedom not found in the physical system. To illustrate 
this, we will continue to use the integration of Eqs. (||) 
as the example upon which we base our outline of these 
integration algorithms. 

As mentioned above, the even and odd degrees of free- 
dom in Eqs. (||) can be written x = (r, u) and y = (q, v). 
But when the force law lets a = a(r) be calculated inde- 
pendently of v, the (r, v) pair of variables are not coupled 
to the (q, u) pair, so this second pair can be dropped from 
the numerical algorithm. This leads to the methodology 
of the staggered leapfrog algorithm, which defines the 
variables it evolves on alternating time levels only. A 



staggered leapfrog integration of Eqs. (0) would, for ex- 
ample, evaluate only r at even steps and only v at odd 
steps. It is then customary to renumber the steps so that 
the even steps are integer values, the odd half integer: 

fn+i = ? n +v n+1/2 At 



fn+3/2 



tn+1/2 



tn+1 



At 



(11) 



where the constant At is the difference between two 
consecutive integer (or half-integer) time levels [E5|. 
The initial conditions are specified by giving (r ,^ 1 ^). 
This method is time-symmetric and symplectic and thus 
avoids nonphysical damping p^ , ^3| . Here, with velocity 
independent Newtonian forces, one has a 2-level method. 
[I.e., only the v components are used in updating the Fs, 
and only the r components in updating the iPs.} This 
leapfrog method gives second order accuracy at the same 
computational cost as the first order Euler method. 

When a depends on v as well as r, a method such as 
extrapolation (e.g., v n+1 = (3/2)v n+1 / 2 - (l/2)u n -V 2 ) 
is needed to calculate v at the integer (or even) levels; 
this generally destroys the time-symmetric nature of the 
algorithm and leads to errors in conserved quantities. 

The Crank-Nicholson technique is an iterative integra- 
tion algorithm p6|. For the system given in Eqs. (||), the 
iteration cycle is initialized by setting 

f n+1 =f n , v n+1 =v n . (12) 

Then data at the half integer levels n + 1 is determined 
via 



-■n+l/2 _ I/y«+l , f>n\ 

2 V ' 

jn+l/2 = l^n+l+^n) 



(13) 



Next, the values of the data at level n + 1 are updated 
according to 



v n+1 = v> 



v n+1 ' 2 At 
a n+l ' 2 At 



(14) 

where a n+1 / 2 has been computed based on the data given 
by Eqs. (|l|). The steps given by Eqs. ( [l3|) and Jl4|) are 
then repeated until the relative differences of the values 
of each of the components of r n+1 and v n+1 between ad- 
jacent iteration cycles have each converged to the desired 
level of accuracy. This method requires the specification 
of initial data r° and v°. It has accuracy 0{At 2 ); it 
is time symmetric if the iteration converges to machine 
precision. 

Both the staggered leapfrog and the Crank-Nicholson 
algorithms avoid the instability arising from an aug- 
mented phase space. One sees in each case that there 
is a single physical pair of values at each timestep. These 
are (f n ,v n+1 / 2 ) for staggered leapfrog and (r n ,v n ) for 
Crank-Nicholson. Other quantities such as v n+1 in 
staggered leapfrog or the half integer values in Crank- 
Nicholson are temporary variables computed from the 
physical data. Thus no extraneous degrees of freedom, 
like the w quantities in 3-level leapfrog, appear. 
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III. DELOUSED LEAPFROG 

The instability in the 3-level leapfrog scheme is ide- 
ally eliminated with a technique that retains as much of 
the symplectic character of the 3-level scheme as pos- 
sible. Aoyagi and Abe (l^j2^] give such a prescription 
to remove the parasitic modes before they grow large 
enough to destroy the calculation. This method relies 
on the Runge-Kutta (RK) method M] and is called the 
RK smoother. (The usual RK algorithm is fourth order; 
however, Aoyagi and Abe do not state what order RK 
scheme is used in their smoother.) They compare it to a 
less effective second order smoother suggested by a col- 
league. We have extended the work of Aoyagi and Abe 
to produce a second order algorithm which requires less 
storage and also allows the use of adaptive (i.e., not con- 
stant in time) timesteps At (see also Q). We term the 
3-level method with this improved smoother "deloused 
leapfrog," since it removes the parasitic modes from the 
calculation. In this section, we give the algorithm for this 
deloused leapfrog method. Ref. [^IJ discusses its proper- 
ties in Hamiltonian problems with applications to highly 
noncircular and highly relativistic geodesies. Those ap- 
plications use the adaptive timestep our delousing routine 
allows. 

The prescription for a delousing step to remove the 
parasitic modes uses a second order Runga-Kutta (RK2) 
algorithm. (The usual RK algorithm is fourth order.) 
Our algorithm proceeds as follows. Referring to the 3- 
level discretization in Eq. (^), we assume that data is 
available on (say) an odd level n and an even level n — 1 . 
First, use RK2 to evolve the data a half-step backward 
from the odd level n and a half-step forward from the 
even level n — 1: 

Step (1) use RK2 with St = -§At 

, . -,n— 1/2 ->7l— 1/2 c ->n ->n 

t0 & etr odd > v odd homr n ,v n 



Step (2) use RK2 with St = + \At 

to get f even 2 , V even 2 from r n ~ 1 ,V n ~ l , 

where, for example, ? n od ^ 2 represents data on level n — | 
that comes from the odd level n. We now have data from 
both odd and even solutions at the same time level n — ^ . 
Eqs. (0) are then utilized to obtain the data f n ~ 1 / 2 and 



n-l 



that contain only the physical mode: 



X ,-m-l/2 



t n-l/2, 
odd / 



Step (3) set f™- 1 / 2 = §( 

— 2 \ v even ~t V odd )■ 

With the parasitic mode thus eliminated at level n — | 
(i.e., with the constraint w = enforced), we now use 
RK2 to step the data a half-step forward and backward: 

Step (4) use RK2 with St' 

to get r n , 

Step (5) use RK2 with St' = -\At' 



-lAt' 



from p»-l/2 ?i p-l/2 



Since RK2 does not suffer from the leapfrog instabil- 
ity (for the same reason as that given in the proceed- 
ing section for the stability of the staggered leapfrog and 
Crank-Nicholson methods), no parasitic mode is intro- 
duced by its use. With this "deloused" data at levels n 
and n — 1, which contain only the physical modes, the 
3-level leapfrog integration is resumed. This prescription 
can be coded so that no more than three levels of data for 
r and v need be kept in memory at any one time, which 
is an important consideration when it is applied to 3 + 1 
dimensional partial differential equations. 

Steps (4) and (5) of this procedure can also be used 
to start the integration. With initial data specified (e.g., 
at t — 0) as values for r 1 / 2 ,?/ 1 / 2 , this procedure pro- 
vides r°,v° and r to start the 3-level algorithm of 
Eqs. (§). 

Notice that the timesteps ±^At' used in Steps (4) and 
(5) to bring the deloused data from level n — \ to lev- 
els n and n — 1 are not required to be the same as the 
timesteps ±^Ai used in Steps (1) and (2). This is an 



improvement on the method of Aoyagi and Abe [|18 22 
and allows the system timestep to be changed whenever 
delousing is performed, with the integration restarted in 
a time symmetric manner. In this paper, all our runs 
were carried out using constant timesteps; see, however, 
Ref. H] for an example of adaptive timesteps applied to 
highly eccentric particle orbits. 

The deloused leapfrog method is not strictly symplectic 
because the delousing steps are taken using RK2, which 
is not a symplectic method, and the suppression Step (3) 
does not preserve volume in the augmented phase space. 
Since delousing steps can in principle be rare, and ideally 
do not disturb the physical solution (when the C(w 2 ) 
terms in Eqs. (|Io| ) are beyond machine precision), the 
failure here of the time symmetric and symplectic prop- 
erties possessed by Eqs. (|j) could have negligible impact. 

To determine when a delousing step is needed, the 
growth of the parasitic mode is monitored. Since Eq. (0) 
shows that this parasitic mode changes sign on each time 
level n, we look for this alternate timestep oscillation in 
some quantity that characterizes the system. The detec- 
tion of such an oscillation in the quantity monitored trig- 
gers a call to the delousing module, which then removes 
the parasitic mode. The goal is to monitor a quantity 
that manifests the instability at detectable levels early 
enough, so that it can be eliminated before it dominates 
the integration. The best quantity to use for the de- 
lousing trigger varies from problem to problem and also 
depends on what level of parasitic mode one is willing to 
tole rate in the solution (see further discussion in Sees. VI 
and[vT|). 

Figure |^ shows the results of our deloused leapfrog in- 
tegration of the same particle equations of motion as 
in Fig. |l|. The delousing trigger chosen in this case 
was an oscillation in a quantity based upon the action 
dJ = \ {pdq — qdp) of Hamiltonian mechanics: 



to get r 



n—l ^tn — 1 



frc 



i-l/2 ^n-1/2 



AJ=|(p' 



P 



L ) • (f n -f n - 1 ) 
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-n-l\ 



■P 



P 



(15) 



where r and p are the particle's position and momentum 
vectors, respectively. The use of this trigger initiated a 
delousing step about once every 405 timesteps. The orbit 
was stable for the duration of our simulation, about ten 
thousand orbits. 



IV. THE 3 + 1 FORMALISM OF NUMERICAL 
RELATIVITY 

The development of the deloused leapfrog algorithm 
described in the proceeding section was motivated by our 
search for a stable and efficient time integration technique 
that could be utilized in CPU and memory intensive, gen- 
eral relativistic simulations of the orbital stability of bi- 
nary neutron stars. This section briefly summarizes the 
3 + 1 form of the Einstein equations upon which such 
numerical relativity simulations are based. Here we con- 
sider only vacuum spacetimes; however, our results apply 
to models with sources as well. 

The Einstein field equations of general relativity are 



G L 



8ttT u 



where G M „ and T^ v 



are the Einstein and 



stress energy tensors, respectively |27j. We set = 
(as is the case for vacuum spacetimes), G = c = 1 and 
our convention is such that Greek indices run from to 
3, Latin indices run from 1 to 3, and repeated raised 
and lowered indices are summed over. The numerical 
solution of the Einstein equations is facilitated by the 

|2? 



ADM decomposition 27 



which divides four dimen- 
sional spacetime into a series of three dimensional, space- 
like slices that are connected by one dimensional, timelike 
curves. This "3+1" split transforms the Einstein equa- 
tions into two sets: the evolution equations, which govern 
temporal changes in the gravitational field variables, and 
the constraint equations, which the field variables must 
satisfy on each spacelike slice. 

The general form of the metric in the ADM formalism 



-{adtf + g ij (dx i + (3 l dt)(dx j + p J dt). (16) 



The vacuum evolution equations for the three-metric g^ 
and the extrinsic curvature Kij are (see, e.g., Q and 
references therein) 



h t = -2aKij + Di[3j + DjPi 



(17) 



and 



K iJtt = -DiDja + ftDtKij + KuDj/3 1 + K tj D^ 1 



+a\R l 



2K tl K + KK, 



(18) 



where t is coordinate time, commas denote partial deriva- 
tives, K = K\, and Di and Rij are the covariant deriva- 
tive and Ricci tensor, respectively, formed from the three- 
metric . The lapse a and the shift vector (3 1 are freely 



specifiable and define the coordinate conditions under 
which g^ and Kij are evolved. The right hand side of 
Eq. (|l7]) may be a function of gij, as dependence on g^ 
arises through the covariant derivative operation when (3 l 
is nonzero or when a is chosen to depend on g^; more- 
over, the right hand side of Eq. ([t8]) is always a function 
of K^. Thus "velocity dependent forces" are generally 
present in this system of equations and therefore 3-level 
leapfrog integrations of them will suffer from the insta- 
blity described in Sec. II. 

The set of constraint equations is comprised of the 
Hamiltonian constraint: 

R- KijK ij +K 2 = 0, (19) 

where R = R\, and momentum constraints: 

Dj{K i3 — g lj K) — 0. (20) 

When the field variables and coordinate conditions for 
the spacetime are functions of time only, the Einstein 
equations [Eqs. ( |l7| ) and (|l8|)] governing the evolution 
of vacuum spacetimes become ordinary differential equa- 
tions: 



dgij 
dt 



-2aKi 



dK„ 



iU a{KKij - 2K a K\). 



(21) 



(22) 



This simplified set of equat ions is the basis of the toy 
codes discussed in Sec. VIA. 



V. GOWDY T A SPACETIMES 

In this section, we introduce the Gowdy T 3 spacetimes 
and then describe three classes of these spacetimes, for 
which exact evolution solutions are known and which we 
have used as relativistic testbeds for the deloused leapfrog 
technique (see Sec. VI.). 

The Gowdy T 3 spacetimes are solutions of the vacuum 
Einstein equations [Eqs. ( |T7| ) and (|l8|)] on the 3-torus, 
in which plane gravitational waves are contained within 
an expanding universe ]30| , Solutions of Gowdy's equa- 
tions have been used for studying the nature of the initial 
cosmological singularity (|Jl^l, and, as is the case in this 
work, for testing numerical codes for solving Einstein's 
equations |3~l) . 

The Gowdy T 3 metric can be written: 



ds 2 = e^l 2 {~e 2T dT 2 + dz 2 ) + e T dw 2 , 



where 



dw z — e r (dx + Q dy) 2 + e p dy 2 



(23) 



(24) 



This form that appears in |T(| is there attributed to Mon- 
crief |32| , whose exploitation of the harmonic map |33 34 



G 



character of the associated Einstein equations could sug- 
gest many equivalent forms. The Gowdy coordinates r, 
A, a, 5, and 9 from [[[o| have here been written as — r, 
—A, x, y, and z. The change in the signs of the time 
coordinate t and of A signify that time increases as the 
universe expands; in most earlier work, the emphasis on 
the study of the initial singularity dictated a time co- 
ordinate that increased as the universe contracted. The 
metric parameters A, P, and Q are functions of z and r 
only and are periodic in z. 

We found the use of an alternate time variable 



t = e T 



(25) 



convenient in that it here makes the (coordinate) veloc- 
ity of propagation of gravitational waves constant. This 
means that with a fixed spatial discretization size Az 
the Courant condition will call for a fixed timestep At 
and thus makes the evolution easier to follow numerically. 
With this time coordinate, the Gowdy metric [Eq. (23)] 
becomes 



ds 



t^ 2 e x / 2 (-dt 2 + dz 2 )+tdw 2 



(26) 



With an ingenious parameterization of the metric sim- 
ilar to Eq. (|24|), Gowdy found that the vacuum Einstein 
evolution equations could be written in terms of P and 
Q alone. With the metric in the form of Eq. (^6|), these 
become 



Pm + * _1 P t - P zz - e 2P {Q 2 t - Q 2 Z )= 
Q,tt + t^Qt - Q,zz + 2(P, t Q,t - P,zQ,z)= 0, 

with the constraint equations specifying A: 
X, z = 2t(P z P, t + e 2P Q, z Q, t ) 



A t = t[P 2 + P 2 



„2P 



(Q 2 t 



,)]• 



(27) 
(28) 

(29) 
(30) 



A. Kasner universe in power-law form 

When the metric parameters P and Q are set to zero, 
the Gowdy metric [Eq. (|26|)1 reduces to the simple diag- 
onal form 

ds 2 = t- 1/2 e x/2 {-dt 2 + dz 2 ) + t{dx 2 + dy 2 ). (31) 

This form of the metric represents a homogeneous, 
anistropically expanding universe with no gravitational 
waves and is equivalent to the form of an axisymmet- 
ric Kasner metric |35| with a different time coordinate. 
A comparison of Eqs. ( |l6| ) and ( |3l"l ) yields the following 
analytic solution for the time evolution of the diagonal 



components of the field variables {gij 
and coordinate conditions: 



for i 7^ j) 



911 = 922 = t , g 33 
Ku = K 2 2 = -jt 1/A , -n-33 

and a = .A/33 = £~ 1/4 



-1/2 



K33 = \t-"' A 



(32) 
(33) 
(34) 

(35) 



where we have used Eq. ( [l7|) to determine Kij in terms 
of gi^ti P l j and a, and have set A, a constant in this case, 
to zero. 



B. Kasner universe in exponential form 

We can write the metric for the same Kasner universe 
with a different time coordinate in an exponential form 



ds 2 = 



l di 2 



5 2t / 3 dx 2 



Q 2t/3 



dy 



2 , e -t 



'Hz 2 



which is essentially the Gowdy metric given in Eq. (23) 
with P = Q = A = 0. A comparison of Eqs. (Bq) and (18) 
gives the following analytic solution for gu (t jand K^Jt) 
(gij = K.^ = for i =/= j) and the coordinate conditions: 



gu = 922 



„2t/3 



533 



-t/3 



-<M1 = ^22 = 



#33 = \e 



-5t/6 



ff = 



and 



= V5 = y/9ii ■ 922 ■ 933 = e 



t/2 



(37) 
(38) 
(39) 
(40) 



where g = det(gry) and we have once again used Eq. ( |T^ ) 
to determine K+j. 

C. Polarized waves in an expanding universe 

The Gowdy metric [Eq. (|26|)1 with Q set to zero, 

ds 2 = t- 1 ' 2 e x / 2 {-dt 2 + dz 2 ) + t(e p dx 2 + e- p dy 2 ), 

(41) 

represents the spacetime of an expanding universe con- 
taining polarized gravitational waves propagating in the 
z-direction. With this metric, the evolution equations 
[Eqs. ( p7| ) and (p8|)] reduce to a single linear equation for 
P: 



Pm + t- x P,t - P,zz = 0. 
The constraint equations become 
\ z = 2tP z P, t 

and 

\, t = t(P 2 + P 2 z ). 



(42) 
(43) 
(44) 



The general solution to Eq. (^2|) is a sum of terms 
of the form Zo(2,Trnt) cos(27rnz) and Zo(2Trnt) sin(27roz), 
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where n is an integer (assuming periodicity of 1 in z) and 
Zq is a linear combination of the Bessel functions Jo and 
Y . We have chosen to study the spacetime based on the 
particular solution 



P = J (27rf)cos(27rz). 



(45) 



A comparison of the metrics given in Eqs. ([l6|) and j4l| ) 
then yields the following exact solution for the time evo- 
lution of the diagonal components of the field variables 
(c/ij = Kij = for i j) and coordinate conditions: 



5ii 



te P , 322 



te P , 333 



t -l/2 e A/2 



K22 
K33 



and 



_l i l/4 e -A/4 e P (1 + t p )t)> 
_lfl/4 e -A/4 e -P (1 _ t p t)) 

£t-vv/4 (t -i_ Ajt)j 

= 0, 



t -l/4 e A/4 



(46) 



(47) 



(48) 



(49) 



where we have again used Eq. ( [L7| ) to determine ify. 
This set of equations is completed by an expression for 
A, which can be derived by using Eq. (Ejj) in conjunction 
with Eqs. © and (||): 



'2TrtJ a (2TTt)J 1 (2Trt) cos 2 (2ttz) 
f 2 7 r 2 t 2 [J 2 (27ri) + J?(2?ri)] 

- H( 2 ^) 2 [Jo(2tt) + Ji 2 (2tt)] - 27rJ (27r)Ji(27r)} 



(50) 



A. Toy code simulations 

The vacuum Einstein equations for spatially homoge- 
neous metrics reduce to the set of ordinary differential 
equations Eqs. and (p^Y. These are similar in form 
to the set of equations Eqs. (^), with the three- metric g^j 
replacing r and the extrinsic curvature Kij replacing v. 
Note that, as in the case of the full Einstein equations 
[Eqs. © and ©], the system of Eqs. © and Eqs. 
(E^) generally contains "velocity-dependent forces," and 
thus 3-level leapfrog integrations of these equations will 
be inherently unstable. We have constructed toy codes to 
solve these equations using each of the integration meth- 
ods discussed in Sec. [n]. The 3-level leapfrog toy code 
is based on the discretization in Eqs. (Q), with the de- 
lousing module based on the Steps (1) - (5) outlined in 
Sec. [II. The staggered leapfrog toy code is based on 
the discretization in Eqs. (pi]); in this method we use 

the extrapolation K™ +1 = (3/2)^" +1/2 - (1/2)X™~ 1/2 
to obtain Kij on the full integer levels p7y . Although 
this extrapolation is accurate to second-order in At, it 
is not time symmetric. Finally, the Crank-Nicholson toy 
code is based on the discretization in Eqs. ([ijj) - (14). 

A toy code simulation of the Kasner universe in power- 
law form from Sec. VA was the initial testbed numeri- 
cal relativity problem to which we applied the deloused 
leapfrog technique. To demonstrate the need for the de- 
lousing modification of the standard 3-level leapfrog tech- 
nique, we first carried out a run with the standard 3-level 
leapfrog technique itself. Eqs. (|3^)-(|3^) provide both the 
initial conditions for this simulation (begun at i = 1 and 
run with a timestep At = 0.1) and the means to measure 
its accuracy via a comparison between the analytically 
and numerically determined values of the field variables. 
Specifically, we use 



VI. NUMERICAL RELATIVITY SIMULATIONS 

We have used the classes of Gowdy spacetimes dis- 
cussed in Sec. [v| as testbeds for numerical relativity sim- 
ulations using the deloused leapfrog, staggered leapfrog, 
and Crank-Nicholson time integration techniques dis- 
cussed in Sec. We have carried out these simulations 
using two types of codes, toy codes (which solve the Ein- 
stein equations Eqs. (|2l|) and j2^ ) for spatially homo- 
geneous spacetimes) and the ADM code developed by 
the Binary Black Hole (BBH) Grand Challenge Alliance 
p6jp] (which solves the full vacuum Einstein equations 
Eqs. p7|) and ( |l8| ) in three spatial dimensions). In this 
section, we present the results of these runs plus efficiency 
analyses that compare the numerical cost effectiveness of 
the various time integration schemes. 



3 

[B 

i=l 



1/2 



(51) 



as a measure of a simulation's accuracy. Here the analytic 
values of the diagonal components gu [given in this case 
by Eq. (|32|)1 are denoted as gf™. Another quantity useful 
in determining the quality of the integration is the size 
of the normalized residual of the vacuum Hamiltonian 
constraint [Eq. (|l9|) 



\R + K 2 



KjjKV 



\R\ +K K] 



(52) 



R\ = R = for the spatially homogeneous Gowdy space- 
times. The results of the 3-level leapfrog simulation are 
presented in Fig. |3|. Panels a and b of this figure dis- 
play both the analytical (solid lines) and numerical (dots) 
solutions for the field variable components gu(t) and 
Kn(t), respectively; panels c and d show the accuracy 
measures H norm {t) and e g (t), respectively. The separate 
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even and odd timestep numerical solutions, character- 
istic of the leapfrog instability described in Sec. 0(cf. 
Eq. (@)), are clearly visible in all four panels of Fig. |3| as 
two dotted branches representing the numerical solution. 
The fact that these two dotted branches are indeed al- 
ternate timestep oscillations of the numerical solution is 
evident in the inset of Fig. |^a, which is an enlargement of 
the numerical solution of 511 for 37.9 < t < 39.2. In this 
inset every value of g±\ has been plotted with an "x" and 
consecutive values have been connected by dashed lines. 

A toy code integration starting with the same initial 
conditions as the 3-level leapfrog run depicted in Fig. || 
was also performed with the deloused leapfrog technique. 
The trigger chosen to initiate the delousing steps in this 
run was a change in the sign of the slope of H norm (t); 
such a sign change is indicative of the alternate timestep 
oscillations discussed in the preceeding paragraph. The 
results of this integration are given in Fig. m Panels a 
and b of this figure display the analytical (solid lines) 
and numerical ("x"s) solutions for g\\(t) and K\\{t), re- 
spectively; the solid lines in panels c and d represent 
\og{H norm {t)) and e 9 (t), respectively. The regular be- 
havior of the evolved quantities shown in Fig. H demon- 
strates that the delousing steps successfully removed the 
parasitic mode from the solutions for the field variables, 
producing an evolution that was stable for the duration 
of the integration. We ran the deloused code a factor of 
25 times longer than the duration of the catastrophically 
unstable 3-level leapfrog simulation. 

Note that in general it is possible to evolve components 
of gij or Kij using the constraint equations [Eqs. ( [l9f - 
pfj| )1 in place of one or more of the evolution equations 
[Eqs. (|T7)-p~8|)] . For integrations of spatially homogeneous 
spacetimes, only the Hamiltonian constraint is mean- 
ingful in this context. In order to investigate the sta- 
bility of such constrained evolutions, we performed 3- 
level leapfrog toy code simulations of the Kasner universe 
in exponential form in which we replaced the evolution 
equation for K 33 [Eq. (||) with i = j = 3] with Eq. (JToJ) . 
Thus -K33 was calculated in terms of the evolved quan- 
tities if n and K22 by imposing the Hamiltonian con- 
straint. These constrained runs still suffered from the in- 
stability under discussion. However, the times at which 
the simulations became catastrophically unstable were 
almost two and a half times longer than in the corre- 
sponding unconstrained runs. 

We also carried out integrations of this model using 
the staggered leapfrog and Crank-Nicholson techniques. 
These runs were all stable, as expected. 



Efficiency analysis of toy code runs 

The efficiency of a stable integration technique is also 
an important factor to consider in evaluating numerical 
methods. Here, we consider the efficiency of a technique 
to be the accuracy level it maintains for a particular nu- 



merical cost. Since the evaluation of the right-hand sides 
of the discretized equations [e.g., Eqs. (H), and jl^ ) 
for the set of equations Eqs. (Q)] is generally the most 
expensive operation in terms of CPU time, we define the 
cost of an integration to be the number of times the right- 
hand sides are computed. 

The results of our efficiency comparison for the toy 
code simulations of the Kasner universe in power-law 
form are displayed in Fig. |5[ For this comparison, we 
ran simulations with each of the three stable integra- 
tion methods; in these simulations the initial conditions 
and evolution duration were identical but the constant 
timestep used during the simulations At was varied from 
run to run. We used the value of e g at the end of the sim- 
ulation as the accuracy measure. The left panel of Fig. || 
gives the final values of log(e g ) as a function of log(Ai) 
and demonstrates that all three techniques are second- 
order accurate (i.e., the slope of log(e g ) versus log(Ai) 
for each method is ~ 2). The right panel of Fig. g gives 
the final values of log(e 9 ) as a function of the numerical 
cost (measured by the number of times both dg^j dt and 
dKij/dt are computed) and provides the most informa- 
tive picture of the efficiency of the different techniques. 
The average number of iterations per timestep for the 
Crank-Nicholson runs ranged between two and three [for 
a convergence criterion of 1.0 x 10 -8 ; see Sec. II]. The 
total number of delousing steps performed during the de- 
loused leapfrog runs ranged from 153 to 181. Note that 
each delousing step adds eight calls to the cost of the 
integration since it requires four calls to the RK2 rou- 
tine (see Sec. Ill), which in turn computes dgij/dt and 
dKij/dt twice. Panel b shows that, for simulations of this 
simple, spatially invariant spacetime, Crank- Nicholson is 
the most efficient of the three integrators. We have used 
least-squares analysis to fit the best straight lines to the 
data points shown in Fig. ||d. This analysis led to the 
following relationship for e g (cost): 



e„ = 10" s cost r ' 



(53) 



The values of the parameters b g and m g for these fits are 
given for each integrator in Table 1. 

We have also done an efficiency comparison of these 
three integration techniques for toy code simulations of 
the Kasner universe in exponential form; the results are 
shown in Fig. ^. The outcome of the efficiency tests for 
these simulations is quite different from that based on 
the Kasner universe in power-law form of Fig. |^. Fig. [3|d 
shows that for relatively low to moderate cost and ac- 
curacy demands, the staggered leapfrog method is the 
most cost effective technique in this case; however, the 
deloused leapfrog method is more efficient when high 
accuracy levels are required. The higher average num- 
ber of iterations per timestep required by these Crank- 
Nicholson runs, which ranged from three for At = 0.0016 
to seven for At = 0.1, may account for the reversal in 
its relative cost effectiveness from the simulations of the 
power-law form (Fig. ||). The number of delousing steps 







taken during the deloused leapfrog runs ranged from 1285 
for At = 0.0016 to 422 for At = 0.1. Thus the deloused 
leapfrog method had to work harder to maintain stable 
integrations of this universe in exponential form than in 
power-law form. We have again used least squares analy- 
sis to fit straight lines to the data in Fig. and produce 
a relation of the form of Eq. (53); this relation is pa- 
rameterized by the values of b g and m g given in Table 
2. 



B. ADM code simulations 

Our study of the deloused leapfrog method stems from 
our search for an efficient technique capable of performing 
numerically stable simulations of the orbital dynamics of 
binary neutron stars. Because such simulations require 
the solution of the full Einstein equations, we wanted to 
test the deloused leapfrog integrator in conjunction with 
the code we plan to use to do these simulations, the ADM 
code developed by the BBH Alliance |3^,|l| . This second- 
order accurate code currently solves the vacuum Einstein 
equations on a Cartesian grid and provides the user the 
choice of utilizing either the standard 3-level leapfrog or 
Crank-Nicholson integration techniques. We have added 
the capability of using the deloused leapfrog integrator 
to the BBH Alliance's ADM code and have used it to 
perform simulations of a Kasner (homogeneous Gowdy) 
expanding spacetime with the power-law coordinate con- 
dition of Sec. V. A and of the expanding Gowdy spacetime 
with polarized gravitational waves of Sec. V.C. Of course 
this code, like the toy code, was ignorant of Gowdy's in- 
genuity which, through parameterizations like gn — te p , 
can reduce some of the Einstein equations to linear equa- 
tions. The Einstein equations are coded in terms of the 
gij and Kij as shown in Eqs. ( |l7| ) and (|l^); the chosen 
coordinate conditions were /3 Z = and a — -^533. They 
involve not only the polynomial nonlincarities manifest 
in these equations, but also the nonlinearities implied 
through a and through the inverse metric when indices 
are raised or covariant derivatives or curvatures are com- 
puted. 

The preliminary testbed used in our ADM code runs 
was a simulation identical to the toy code runs of the 
Kasner power-law metric. The development of the insta- 
bility in the ADM code's 3-level leapfrog run replicated 
its development in the toy code run. The ADM code's 
deloused leapfrog run successfully removed this instabil- 
ity in the same manner as in the toy code run, with the 
delousing steps triggered at the same temporal intervals 
in both simulations. 

To test the behavior of the deloused leapfrog method 
for partial differential equations with spatially varying 
terms, we carried out simulations of the polarized Gowdy 
spacetime of Sec. V.B. Equations (45)-(p0[) yield both 
initial conditions for simulations of this spacetime and 
exact solutions with which to compare the results of such 



simulations. 

Our ADM code polarized Gowdy simulations began 
at t = 1 and were run with periodic boundary condi- 
tions over the interval — | < z < i and a grid spac- 
ing Az = 1/62. Because the vacuum Einstein equa- 
tions are partial differential equations, the size of the 
timestep that can be taken in the integration is restricted 
by the Courant condition pqj , which ensures that infor- 
mation cannot propagate across more than a single grid 
zone in one timestep. For the polarized Gowdy met- 
ric of Eq. (|4l]), this condition is equivalent to enforcing 
At = C Az, where Az is the (uniform) grid spacing and 
the Courant factor C < 1. In the runs presented here, 
we chose C = 0.3. The initial ADM code simulation 
was performed with the 3-level leapfrog integrator and, 
as expected, was unstable. The results of the integra- 
tion are presented in Figs. |?] and |[ The evolution of 
the value of gn at the center of the grid is shown in 
panel a of Fig. 0; panels b, c, and d of this figure dis- 
play, respectively, H norm (t), e g (t), and ~e~t{t). Here bars 
denote (spatial) averages over the grid. H norm {t) is com- 
puted by dividing the spatial average of the numerator 
of Eq. ([32]) by the spatial average of the denominator of 
Eq. (p2). The additional measure of the accuracy of the 
simulation, et, is defined by 

e t (x,y,t) = t~ l \t- [g 11 (x,y,x)g 2 2(x,y,x)} 1/2 \ . (54) 

The usefulness of et as an error estimate arises because, 
according to the analytic solution of Eq. (46), (711322 = t 2 - 

The separate even and odd timestep solutions, indica- 
tive of the instability, can clearly be seen in the plots 
of H norm (t), e g (t), and et(t); however, the separate so- 
lutions are not yet visible in the plot of gn at t = 5.8. 
As shown in Fig. |^, they do appear in gn at the grid 
center later in the evolution, as the instability begins to 
overwhelm the computation. 

We then evolved the same polarized Gowdy initial data 
with the ADM code using the deloused leapfrog integra- 
tor. Because H norm {t) and e g (t) oscillated in our simu- 
lations of this spacetime, they were not used as a basis 
for the trigger that initiated the delousing steps. [These 
fluctuations are due to the complex oscillatory nature 
of g^ and Kij in this spacetime and are not related to 
the alternate timestep oscillations caused by the insta- 
bility; in fact, such fluctuations are also present in the 
Crank-Nicholson simulations of this spacetime.] Instead, 
because e t behaved monotonically (see Fig. ^|d), a change 
in the sign of its temporal slope was used as the delous- 
ing trigger. As can be seen in Fig. ^, the removal of 
the parasitic mode during the delousing steps taken in 
this simulation eliminated the presence of large alternate 
timestep oscillations and allowed for a stable integration. 



s of ADM code runs 



Efficiency 



We have also used simulations of this polarized Gowdy 
spacetime to evaluate the efficiency of the deloused 
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leapfrog algorithm. However, in this case its performance 
could only be compared with that of the Crank-Nicholson 
technique, as the staggered leapfrog method has not been 
implemented in the ADM code. (The reason for this is 
that the memory requirements of the staggered leapfrog 
method would exceed those of the other two methods if 
adaptive mesh refinement were to be used.) In this ef- 
ficiency comparison, Az and At are reduced in tandem 
from run to run, with C held constant at 0.3. The value of 
e t at the end of the simulations was used as the measure 
of accuracy upon which to base the efficiency analysis in 
this case. 

The results of this analysis are displayed in Fig. [!(]. 
The left panel of this figure gives the final value of et 
as a function of the grid spacing Az, which was varied 
from 1/126 to 1/62 to 1/30. Straight lines fit throught 
these data points have slopes ~ 2; this demonstrates 
the second-order accuracy of the deloused leapfrog and 
Crank-Nicholson methods. The right panel contains a 
plot of e t versus cost and indicates that the deloused 
leapfrog integrations of this spacetime were about five to 
eight times more efficient than those carried out with 
the Crank-Nicholson technique. The average number 
of iterations per timestep for the Crank-Nicholson runs 
ranged from five for Az = 1/126 to eight for Az — 1/30. 
The total number of delousing steps taken was relatively 
constant in the three deloused leapfrog runs (six for 
Az = 1/126; five for Az = 1/62; and six for Az = 1/30). 
The least squares straight line fits to the data in Fig. nOp 
can be transformed, in this case, to relations of the form 



e t = 10 ht cost 



(55) 



the parameters bt and m t determined by these fits are 
given in Table 3. 



VII. CONCLUSIONS 

The purpose of this paper is to alert the community 
to the existence of the instability inherent in standard 3- 
level leapfrog integrations of Einstein's equations and to 
demonstrate that the proposed delousing modification to 
the standard algorithm can efficiently cure this instabil- 
ity. To this end, we have used three classes of testbed 
solutions. 



In Sec. Ill 



we show calculations of highly 
relativistic circular geodesic orbits calculated in Carte- 
sian coordinates. 



In Sec. VI 



we show evolutions of ho- 
mogeneous expanding cosmologies and polarized gravita- 
tional waves in an expanding Gowdy spacetime. In all of 
these cases, the deloused leapfrog algorithm removed par- 
asitic modes from the numerical solution of the Einstein 
equations that were excited to instability in traditional 
3-level leapfrog simulations of these spacetimes, and thus 
allowed for their stable evolution. 

The numerical efficiency (i.e., the accuracy level main- 
tained for a particular numerical cost) of the deloused 
leapfrog integrator was compared to the efficiencies of 



two other stable integration methods, staggered leapfrog 
and Crank-Nicholson. We have defined numerical cost as 
the number of times the right hand sides of both Ein- 
stein evolution equations (i.e., g+jj and Kij,t) are com- 
puted during a simulation. Thus cost in this case is a 
measure of a simulation's CPU expense. Note that all of 
the simulations presented in this paper were carried out 
with constant timesteps. 

The first testbed for this efficiency analysis was a toy 
code (which solves the spatially invariant Einstein equa- 
tions [Eqs. ( pi) ) and ((p2])]) simulation of a spatially ho- 
mogeneous Gowdy spacetime yielding a Kasner universe. 
With the coordinate condition choice a — <Jgla this gives 
a power-law analytic solution for the metric components. 
For this simple problem, Crank-Nicholson was the most 
efficient of the integrators (see Fig. [|). The results were 
different, however, when the testbed was changed to a toy 
code simulation of the same Kasner universe with a dif- 
ferent time coordinate choice a = ^fg = \Jg~\\ • g22 • 333, 
yielding in the analytic solution an exponential form for 
the metric components. In that case, the staggered and 
deloused leapfrog techniques were more cost effective (see 
Fig. ||), as the rapid evolution of the spacetime caused the 
iterative Crank-Nicholson technique to require a larger 
number of iterations per timestep. 

The final testbed used in our efficiency analysis was 
a simulation of an expanding Gowdy spacetime con- 
taining polarized gravitational waves. These simula- 
tions required the solution of the complete vacuum Ein- 
stein equations and were carried out with the BBH Al- 
liance's ADM code pqjjfl, in which the standard (un- 
stable) 3-level leapfrog and the (stable) Crank-Nicholson 
integration methods had previously been implemented. 
We modified this code to allow the use of the deloused 
leapfrog scheme. Because the staggered leapfrog method 
has not been implemented in the ADM code, the efficien- 
cies of only the Crank-Nicholson and deloused leapfrog 
integrators were evaluated in this case. The deloused 
leapfrog integrations of this spacetime were five to eight 
times more cost effective than the Crank-Nicholson runs 
(see Fig. 0). 

Thus, the results of our testbed simulations indicate 
that the deloused leapfrog algorithm is an effective and 
efficient cure for the 3-level leapfrog instability. Further 
evaluation of this algorithm, via its use in more com- 
plex problems in numerical relativity, such as a contract- 
ing Gowdy universe containing unpolarized gravitational 
waves, and other fields, would serve to confirm the ro- 
bustness of the method and provide insight into its cost 
effectiveness in different numerical scenarios. 

One aspect of the deloused leapfrog algorithm that has 
the potential to alter the conclusions of such efficiency 
analyses is the choice of delousing trigger. Based on our 
experience, we suspect that choice of a trigger which ini- 
tiates an excessive number of delousing steps will degrade 
the accuracy of a simulation to some degree. If this is the 
decrease in the average interval between delous- 
ing steps would not only increase the cost of the run, but 
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would also decrease its numerical accuracy somewhat. 
On the other hand, an increase in the delousing interval 
would allow the parasitic mode in the numerical solution 
to grow to higher levels. For example, had a change in the 
sign of the temporal slope of gu(t) been used as the de- 
lousing trigger in the deloused leapfrog simulation of the 
polarized Gowdy spacetime, the parasitic mode would 
likely have grown to a greater extent between delousing 
intervals, as its presence became sizeable in gn(t) rather 
late in the evolution (see Figs. and ||). Thus the choice 
of delousing trigger may involve a trade-off between the 
loss of some degree of accuracy introduced into the com- 
putation via the delousing steps and the degree to which 
the parasitic mode is permitted to grow between these 
steps. 

In conclusion, we have demonstrated that the deloused 
leapfrog algorithm permits the stable numerical evolution 
of simple vacuum spacetimes. In addition, our results 
suggest that deloused leapfrog may be a better integra- 
tion technique than Crank-Nicholson to employ in com- 
plex numerical relativity simulations, as this new algo- 
rithm was more cost effective than the Crank-Nicholson 
method in our simulation of a spatially varying space- 
time. 
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FIG. 1. The numerical integration of the geodesic equa- 
tions for a particle in the Schwarzschild spacetime using the 
3-level leapfrog technique. The geodesic equation was solved 
in rectangular coordinates from a 3D Hamiltonian (see J2lj]). 
The particle was given initial conditions such that it should 
remain on a circular orbit of radius ro = 10M and have an 
orbit period of 199Af. In each frame, the data is plotted on 
every twenty-third timestep (At = 0.1M). The instability 
manifests as the solutions on odd (circles) and even (trian- 
gles) timesteps diverge. Although the integrator appears to 
be stable at early times, the parasitic mode is present from 
the beginning, on a much smaller scale than is used in these 
plots. The integrator failed and the code crashed after ~ 6 
orbital periods, (a) The particle orbit in the x-y plane. (6) 
The magnitude of the particle's position vector as a function 
of time. 



FIG. 2. Same as Fig. |l| except that the integration was 
carried out using the deloused leapfrog method. The orbit is 
now stable for the ~ 10, 000-orbit duration of the simulation. 
The timestep was the same as that used in Fig. |l|, but the data 
is plotted only every 2001 timesteps. Note that even though 
only ~ 1 point per orbit (each orbit is aproximately 2000 
timesteps) are shown, there are over 10,000 points plotted in 
the figure. Both odd and even solutions lie directly on top of 
each other, filling out the orbit track. The delousing module 
was applied on average once every 405 timesteps, hence this 
figure was produced at a 2 percent increase in CPU cost per 
orbit over the simulation shown in Fig. [I]. 



FIG. 3. The results of the unstable 3-level leapfrog toy code 
integration of a Kasner universe in power-law form are pre- 
sented here. The numerical (dots, with every other pair of 
even and odd timestep values plotted) and analytical (solid 
lines) solutions for <7n(t) and K\\(t) are given in panels a and 
b, respectively. The inset in panel a is an enlargement of the 
numerical solution for gu(t) in the range 37.9 <t< 39.2, in 
which all data points have been plotted and connected with 
dashed lines to emphasize the large alternate timestep os- 
cillations of this unstable solution. The numerical accuracy 
measures H norm (t) and e g (t) are shown in panels c and d, 
respectively; for the sake of clarity, only every other pair of 
even and odd timestep values has been plotted. 



FIG. 4. This figure depicts the stable deloused leapfrog, 
toy code integration begun with the same initial conditions 
as the unstable 3-level leapfrog simulation presented in Fig. ^. 
The numerical ("x"s, with every 201st point plotted) and an- 
alytical (solid line) solutions for gu(t) and Kn(t) are given 
in panels a and b, respectively. The accuracy measures 
\og(H n0 rm(t)) and e g (t) are shown in the lower panels c and 
d, respectively. 



FIG. 5. The results of the efficiency analysis of toy code 
simulations of a Kasner universe in power-law form, repre- 
sented by the metric given in Eq. fell), carried out with dif- 
ferent integration techniques are shown here. The final values 
of log(e s ) are plotted versus log(Af) in the left panel and ver- 
sus log(cost) in the right panel. Values from Crank-Nicholson 
runs are marked with "x"s, those from deloused leapfrog runs 
are marked with triangles, and those from staggered leapfrog 
runs are marked with squares. The cost is defined as the num- 
ber of times the right-hand sides of the discretized equations 
are evaluated. The slopes of straight lines fit through data in 
(a) are ~ 2, indicating that the numerical techniques are all 
accurate to second-order in At. 



FIG. 6. The same quantities and notation as in Fig. |5| but 
for toy code simulations of a Kasner universe in exponential 
form, described by the metric of Eq. (Bra). 



FIG. 7. The results of the unstable 3-level leapfrog ADM 
code integration of the expanding universe containing polar- 
ized gravitational waves are presented here. The numerical 
(dots, with every fourth pair of even and odd timestep data 
points plotted) and analytical (solid line) solutions for gn(t) 
at the grid center are given in panel a. The numerical ac- 
curacy measures H norrn , e s , and e t are plotted as functions 
of time in panels b, c, and d, respectively; again, only every 
fourth pair of even and odd timestep values has been plotted. 
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FIG. 8. The extended evolution of the metric component 
gn, at the grid center, from the unstable 3-level leapfrog sim- 
ulation presented in Fig. (?) is shown here. The notation is the 
same as that of panel a in Fig. ^ (except that every other pair 
of even and odd timestep values is plotted) , but the duration 
of the evolution has been extended to exhibit the growth of 
the parasitic mode, as evidenced by the appearance of the 
even and odd timestep branches of the numerical solution. 

FIG. 9. This figure depicts the stable deloused leapfrog, 
ADM code integration begun with the same initial conditions 
as the unstable 3-level leapfrog simulations shown in Figs. [?] 
and |^. The numerical ("x"s, with every eleventh point plot- 
ted) and analytical (solid line) solutions for <?n(t) at the cen- 
ter of the grid are given in panel a. Panels b, c, and d present 
the evolutions of the accuracy measures H norrn , e g , and e t , 
respectively. 

FIG. 10. The results of the efficiency analysis of ADM 
code simulations of an expanding universe containing polar- 
ized gravitational waves, represented by the metric of Eq. 
([Il]), are shown here. The final values of log(e t ) are plot- 
ted as a function of the logarithm of the grid spacing Az 
(At = 0.3 Az) in the left panel and as a function of the log- 
arithm of the numerical cost in the right panel. Triangles 
mark the data points from deloused leapfrog runs; "x"s mark 
those from Crank-Nicholson runs. Straight lines fit through 
data points in the left panel have slope ~ 2, indicating that 
the numerical techniques are both accurate to second-order 
in At. 



TABLE I. Parameters for fits to accuracy versus cost data 
in Fig. | 



Method 


b a 


m g 


deloused leapfrog 


6.6 


-2.0 


Crank-Nicholson 


6.8 


-2.2 


staggered leapfrog 


6.7 


-2.0 


TABLE II. Parameters for fits to 


accuracy versus 


cost data 


in Fig. | 






Method 


b 3 


m g 


deloused leapfrog 


11 


-3.1 


Crank-Nicholson 


9.0 


-2.6 


staggered leapfrog 


5.8 


-1.9 


TABLE III. Parameters for fits 


to accuracy versus cost 


data in Fig. [lo] 






Method 


bt 


m t 


deloused leapfrog 


4.8 


-2.1 


Crank-Nicholson 


9.6 


-2.8 
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